library(ggplot2)
require(data.table)
library("dplyr")
library("ggpubr") 
require(tidyverse) #for easy data manipulation and visualization
require(caret) #for easy machine learning workflow
library(ncdf4)
library(mgcv) 

wd1 <- "C:\\Users\\yfan\\...\\NFood_replication\\"
case = c("RCP45","RCP45_85LU","RCP85","RCP85_45CO2","RCP85_SAI","RCP85_MSB","RCP85_CCT")
casename = c("RCP45",expression('RCP45'[(SSP5LU)]),"RCP85",expression('RCP85'[(RCP45CO2)]),"RCP85+SAI","RCP85+MSB","RCP85+CCT")


#=======================
#Prepare annual crop yield data
yield = read.csv(paste0(wd,"/data/yield_data/cropyield_annual.csv"), stringsAsFactors=FALSE, na.strings =c("NA"," ","9.96921e+36"))
yield = as.data.table(yield)
yield = yield[Year!=2100,] #remove year 2100, which has a sharp jump due to termination of SRMs
nrow(yield)/length(pft)/length(case) #94 years
yield[Year<2025 & Case=="RCP85_SAI",] #make sure NA values are read correctly


#save yield data per crop as named object or data.table
#add rainfed/irrigated as a dummy variable
yield6crop=yield8crop = data.table() #6/8 crops whether/not merge Tropical Corn, Tropical Soybean
for (i in 1:8) {
  rainfed = pft[seq(1,15,2)][i]
  irrigated = pft[seq(2,16,2)][i]
  yield0 = yield[PFT==rainfed,]
  yield0$Irrigated = 0 #dummy variable, 0=rainfed, 1=irrigated
  yield1 = yield[PFT==irrigated,]
  yield1$Irrigated = 1 #dummy variable, 0=rainfed, 1=irrigated
  yield.tmp = rbind(yield0, yield1)
  yield.tmp$Crop.name = crops[i] #use same crop name for rainfed and irrigated
  yield.tmp$Crop = yield.tmp$Crop.name #rename column variable
  yield.tmp[,Crop.name:=NULL] #remove column
  yield8crop = rbind(yield8crop, yield.tmp)
  if (i<=6) {yield6crop = rbind(yield6crop, yield.tmp)}
}
unique(yield6crop$Crop)
unique(yield8crop$Crop)
#data.table sees column names as if they are variables
yield6crop[Crop=="Corn",]$Yield..MtC.year. = yield8crop[Crop=="Corn",]$Yield..MtC.year. + yield8crop[Crop=="Tropical Corn",]$Yield..MtC.year. #temperate+tropical corn
yield6crop[Crop=="Soybean",]$Yield..MtC.year. = yield8crop[Crop=="Soybean",]$Yield..MtC.year. + yield8crop[Crop=="Tropical Soybean",]$Yield..MtC.year. #temperate+tropical soybean
yield6crop[Crop=="Corn",]$Area..Mha. = yield8crop[Crop=="Corn",]$Area..Mha. + yield8crop[Crop=="Tropical Corn",]$Area..Mha. #temperate+tropical corn
yield6crop[Crop=="Soybean",]$Area..Mha. = yield8crop[Crop=="Soybean",]$Area..Mha. + yield8crop[Crop=="Tropical Soybean",]$Area..Mha. #temperate+tropical soybean
yield6crop[,PFT:=NULL]
#verify the values are correct
mean(yield8crop[PFT==17 & Case=="RCP45" & Year >=2090 & Year <= 2099]$Yield..MtC.year.) #rainfed temperate corn
mean(yield8crop[Crop=="Corn" & Case=="RCP45" & Year >=2090 & Year <= 2099]$Yield..MtC.year.) #average of rainfed and irrigated temperate corn
mean(yield6crop[Crop=="Corn" & Irrigated==0 & Case=="RCP45" & Year >=2090 & Year <= 2099]$Yield..MtC.year.) #rainfed temperate corn+tropical corn
head(yield6crop[380:400,])
nrow(yield6crop)/6/7/2 #6 crops x 7 cases x 94 years x 2 (irrigated dummy 0/1)

#merge rainfed+irrigated yield for 6 crops
yield.all = yield6crop[Irrigated==0,]
yield.all$Yield..MtC.year. = yield6crop[Irrigated==0,]$Yield..MtC.year. + yield6crop[Irrigated==1,]$Yield..MtC.year.
yield.all$Area..Mha. = yield6crop[Irrigated==0,]$Area..Mha. + yield6crop[Irrigated==1,]$Area..Mha.
yield.all[,Irrigated:=NULL] #remove 'Irrigated'
#yield.all = yield6crop #keep Irrigated/Rainfed
yield.all[Crop=="Corn", Crop:="Maize"]
yield.all[Crop=="Soybean", Crop:="Soy"]
mean(yield.all[Crop=="Maize" & Case=="RCP45" & Year >=2090 & Year <= 2099]$Yield..MtC.year.) #rainfed + irrigated & temperate + tropical corn

saveRDS(yield8crop,  file = paste0(wd1,"/data/Rdata/yield8crop.all.Rds")) #for use in Fig.3.
saveRDS(yield.all,  file = paste0(wd1,"/data/Rdata/yield6crop.all.Rds"))




#=================Global trend of crop yield==============  

#save two color-blind-friendly palettes, one with gray, and one with black
#http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/
# The palette with grey:
# cbPalette <- c("#999999", "#009E73", "#D55E00", "#FF9999", "#E69F00",  "#0072B2", "#56B4E9", "#CC79A7", "#F0E442", "#9999CC", "#66CC99")
# barplot(height=rep(1,length(cbPalette)),col=cbPalette) #display colors
# # The palette with black:
# cbPalette <- c("#000000", "#009E73", "#009E73", "#D55E00", "#D55E00", "#CC79A7",  "#56B4E9", "#F0E442",  "#E69F00", "#0072B2" )
# barplot(height=rep(1,length(cbPalette)),col=cbPalette)

# cbPalette <- c("#26828EFF", "#26828EFF", "#440154FF", "#440154FF", "#CC6600", "#56B4E9", "#B4DE2CFF")
# barplot(height=rep(1,length(cbPalette)),col=cbPalette) #display colors
casecolor2 <- c("#440154FF", "#440154FF", "#7E6148FF", "#7E6148FF", "#CC79A7", "#8491B4FF", "#91D1C2FF")
barplot(height=rep(1,length(casecolor2)),col=casecolor2) #display colors


#==========plot trend of global total yield (ha)
yield.all <- readRDS(paste0(wd1,"/data/Rdata/yield6crop.all.Rds"))

library(forecast)
yield.gl <- copy(yield.all)
yield.gl <- yield.gl[, .(Prod=sum(Yield..MtC.year.), Area=sum(Area..Mha.)), by =.(Case,Year)]
yield.gl <- yield.gl[, Yield:=Prod/Area/0.45*0.85] #convert yield to tonnes/ha 
yield.gl <- yield.gl[, Yield.ma:=ma(Yield, 5), by=.(Case)]
#head(yield.gl[Case=="RCP85"]); head(df.5case.wm[Scenario=="85"]) #comparable global yield processed by NCL and R scripts

cropname = c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")


#======linegraph of global crop yields  (all 7 cases)
case = c("RCP85","RCP85_45CO2","RCP45","RCP45_85LU","RCP85_SAI","RCP85_MSB","RCP85_CCT")
casename = c("RCP85",expression('RCP85'[('45CO2')]),"RCP45",expression('RCP45'[('SSP5LU')]),"SAI","MSB","CCT")
override.linetype <- c("solid", "longdash",  "solid", "longdash", "solid", "solid", "solid")
override.size <- c(0.6, 0.5, 0.6, 0.5, 0.5,0.5,0.5)
#use mutate to internally reorder yield.all$Crop based on specified levs
p0 <-  yield.gl %>% 
  mutate(Case = factor(Case, levels=case)) %>%
  ggplot(aes(Year, Yield.ma, group=Case)) + 
  geom_line(aes(color = Case, linetype=Case, size=Case) )+ #use Case as group and set override values in manual scale
  #geom_point(aes(color = Case, shape=Shape), size = 0.5) +
  #change linetype, color, size mannually
  scale_linetype_manual(values=override.linetype, guide = FALSE)+ 
  scale_size_manual(values=override.size, guide = FALSE) + #turn off legend for size
  scale_colour_manual(name = "Case", labels = casename,  values = casecolor2) + #, alpha(cbPalette[5:8],0.5))) +
  scale_x_continuous(name="",breaks =c(2006, 2025, 2050, 2075, 2099)) +
  scale_y_continuous(expression("tonnes ha"^-1*"")) #, sec.axis = sec_axis(~ . * 1., name = "Area (Mha)"))
#merge the legends for color and size 
p0 <- p0 + guides(colour = guide_legend(title=NULL,  keywidth = 1.5, keyheight=1, label.position = "right", label.hjust = 0,
                                        override.aes = list(size = override.size*0.8, linetype = override.linetype )))
##remove grid and background or remove all
p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) 
p0 + ggtitle(expression('('*bold(f)*') Global crop yield')) + 
  theme(plot.title = element_text(hjust = 0.5,size=15)) + #theme(strip.text.x = element_text(size = 10)) +
  theme(legend.text = element_text(size=11), legend.position = c(0.24, 0.75), legend.background = element_blank()) + #legend.title = element_blank(), 
  theme(axis.text = element_text(size=12), axis.title=element_text(size=13))

ggsave(paste0(wd1,"/figures/main/Fig1f.GloalYield_Trend.svg"), device ="svg", units = "in", width = 4, height = 3.65, dpi = 600)






#================supplement figure S1==================

library("forecast")
library("zoo")
###show 5-year running mean
#yield_df <- copy(yield.all.cor) #use bias corrected data for soy
yield_df <- copy(yield.all) #use original global data
yield_df <- yield_df[, Prod.ma:=ma(Yield..MtC.year.,5)/0.45*0.85, by="Case,Crop"] #Million tonne per year
yield_df <- yield_df[, Yield.ma:=ma(Yield..MtC.year./Area..Mha.,5)/0.45*0.85, by="Case,Crop"] #convert to tonne/ha by carbon:biomass ration and harvest efficiency
yield_df <- yield_df[c("RCP45","RCP85"), Area:=Area..Mha., on="Case"] #show area of SSP2/SSP5, other cases filled with NA

cropname = c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")

#In order to change the order of the panels, you need to change the factor levels for your faceting variable Crop
levs <- c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")
override.linetype <- c("solid", "longdash",  "solid", "longdash", "solid", "solid", "solid")
override.size <- c(0.8, 0.6, 0.8, 0.6, 0.6,0.6,0.6)
#use mutate to internally reorder yield.all$Crop based on specified levs
p0 <-  yield_df %>% 
  mutate(Crop = factor(Crop, levels=levs, labels=c("(a) Maize","(b) Sugarcane","(c) Wheat","(d) Rice","(e) Soy","(f) Cotton"))) %>%
  mutate(Case = factor(Case, levels=case)) %>%
  ggplot(aes(Year, Yield.ma, group=Case)) + 
  geom_line(aes(color = Case, linetype=Case, size=Case) )+ #use Case as group and set override values in manual scale
  #geom_point(aes(color = Case, shape=Shape), size = 0.5) +
  #change linetype, color, size mannually
  scale_linetype_manual(values=override.linetype, guide = FALSE)+ 
  scale_size_manual(values=override.size, guide = FALSE) + #turn off legend for size
  scale_colour_manual(name = "Case", labels = casename,  values = casecolor2) + #, alpha(cbPalette[5:8],0.5))) +
  #for crop area change
  # geom_line(aes(Year, Area, linetype=Case), color="#999999", size = 0.8 ) + 
  # geom_vline(aes(xintercept=tp.year), linetype="dotted", color="#999999", size = 0.8 ) +
  # geom_text(aes(x=tp.year,y=rep(c(400,50,350,400,130,15),each=94*7), label=tp.year, group=Case), color="gray3") +
  #scale_linetype_manual(values=c("solid", "dashed"), guide = FALSE)+ #can only use this scale for one time
  #use free scales for each panel
  facet_wrap(~Crop, nrow=3, scales='free_y') +  
  scale_x_continuous(breaks =c(2006, 2025, 2050, 2075, 2099)) +
  scale_y_continuous(expression("Yield (tonnes ha"^-1*")")) # , sec.axis = sec_axis(~ . * 30, name = "Area (Mha)"))
#merge the legends for color and size 
p0 <- p0 + guides(colour = guide_legend( keywidth = 2, label.position = "right", label.hjust = 0,
                                         override.aes = list(size = override.size, linetype = override.linetype )))
##remove grid and background or remove all
p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) 
#p <- p + xlab("Year") + ylab("Yield change (MtC/year)")
p0 <- p0 + theme(strip.text.x = element_text(size = 13)) +
  theme(legend.text = element_text(size=14), legend.title = element_blank())+ 
  theme(axis.text = element_text(size=13), axis.title=element_text(size=14))
# Change the color palette
#p0 <- set_palette(p0, palette = cbPalette[-1])
p0 +   theme(axis.title.y.left = element_text(color = "#000000"),
             axis.title.y.right = element_text(color = "#999999"),
             axis.text.y.right = element_text(color = "#999999"))

ggsave(paste0(wd1,"/figures/supplementary/Fig.S1.CropYield_Trend.svg"), units = "in", width = 10, height = 8, dpi = 600)





















